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Abstract 

Dynamic simulators model systems evolving over time. Often, it operates iteratively over fixed number of time-steps. 

The output of such simulator can be considered as time series or discrete functional outputs. Metamodeling is an 

^^ eff'ective method to approximate demanding computer codes. Numerous metamodeling techniques are developed for 

^S| simulators with a single output. Standard approach to model a dynamic simulator uses the same method also for 

^ multi-time series outputs: the metamodel is evaluated independently at every time step. This can be computationally 

Qh demanding in case of large number of time steps. In some cases, simulator outputs for diff'erent combinations of input 

parameters have quite similar behaviour. In this paper, we propose an application of shape invariant model approach 

to model dynamic simulators. This model assumes a common pattern shape curve and curve-specific diff'erences in 

amplitude and timing are modelled with linear transformations. We provide an efficient algorithm of transformation 

parameters estimation and subsequent prediction algorithm. The method was tested with a CO2 storage reservoir case 

using an industrial commercial simulator and compared with a standard single step approach. The method provides 

satisfactory predictivity and it does not depend on the number of involved time steps. 
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1. Introduction 

Simulation models are used nowadays in many industrial applications to predict and analyze the behaviour of 
VO complex systems. A simulator is a complex computer code embedding the physical laws governing the physical 

00 system under investigation. The input of such simulators can be adjustable or uncontrollable parameters which are 

only partially known and thus are aff'ected by uncertainty. Uncertainty analysis of numerical experiments is used 
"^ to assess the confidence of the model outcomes which are then used to make decisions 1 1 1. Here, we focus on a 

^^ particular type of dynamic simulators that are used to make predictions in the future. These simulators are based 

typically on finite element/finite diff'erence codes used for instance for the simulation of flows and transfers in porous 
media. Industrial applications using this type of simulators are for instance hydrocarbons reservoir forecasting and 
carbon dioxide (CO2) underground storage ||2|[3l. Such applications involve very complex numerical codes with a 
large number of inputs and with high uncertainty derived from an incomplete knowledge of subsurface formation O. 
;^ The uncertainty on the simulator output is usually assumed to be mostly due to the propagation of the uncertainty on 

the input. Nevertheless, modelling errors can also play a role. 

The simulator models a multi-phase 3-D fluid flow in heterogeneous porous media, operating over fixed number of 
time-steps. The typical output of such simulators consists of a sequence of outputs at diff'erent time-steps. Therefore, 
it represents time series related, for instance, to a recovery rate for a given well or a group of wells. It can be also a 
spatial output such as a pressure map or a saturation map also evolving with time. Here we focus on ID time series 
output which can be typically measured in a well. 

Let us consider the output of a simulator as a deterministic function Y(0 = F(x, t), where x g Q c R^ is a 
configuration of preselected input parameters and < r < T refers to a time step. Y(0 is a time dependent output, e.g. 
oil production rate or reservoir pressure. 

The function F : R^ x [0, T] ^ R models the relationship between the input and the output of the numerical simu- 
lator. Methods such as Monte Carlo ones can be used to propagate uncertainty. However, each model evaluation being 
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generally very time-consuming this approach can be unpractical in industrial applications. Moreover, the probability 
distribution of the input may vary as new data comes in or within a process of dynamic data assimilation. Therefore, a 
response surface approach (also referred to a metamodel approach) can be more useful in such applications |[3]|5l[6]|. 
The advantage of a metamodel is that is fast to evaluate and it is designed to approximate the complex computer 
code based on a limited number of simulator evaluations. These evaluations of the simulator are taken at some well 
chosen input configurations also called the training set or the experimental design. Numerous experimental designs 
have been proposed and many of them are quite sophisticated. Here, we use Latin Hypercube designs Q for their 
good space filling properties. Usually, it can be coupled with other criteria such as the maximin design (maximum 
minimum distance) 1 8 , 9|. In this work, we focus on Gaussian process (GP) based metamodels also known as kriging 

uniiiiisi. 

The aim of this work is to propose a new method to address multi-time series outputs. For such dynamic simulator 
a standard approach assumes a single step GP model for each time step. This basic approach can be computationally 
intensive for a high number of time steps and when the construction of a single GP model is time consuming (i.e. when 
the size of the training set and the number of variables is large). Therefore, the procedure can become unpractical in 
some industrial applications. 

The problem of metamodeling for dynamic simulators was recently investigated by diff'erent authors and several 
principal approaches can be distinguished. A first possible approach is GP metamodelling considering the time steps 
as the model additional input parameter |[3l|TT][T2]|. This approach is easy to implement, however if we want to take 
into account all the information from an experimental design at every time steps, the size of new experimental design 
is multiplied by the number of time steps. It results to matrices of high dimensions that can lead to some numerical 
problems in case of large size of original experimental design or in case of high density of the steps in the time scale. 
Conti et al. (2009) |[T3]| developed an iterative approach to build a model of dynamic computer code, assuming that the 
model output at a given time step depends only on the output at the previous time step. To reduce the dimensionality of 
the problem, we can also represent the given observations as truncated functions by some selected basis functions with 
the following GP modelling for the coefl&cient of the selected representation. Bayarri et al. (2007) |[T4l introduced 
wavelet decomposition. Campbell et al. (2006) |T5l, Higdon et al. (2008) 1 16|| and Lamboni et al. (2009) |17| 
suggested application of principal component decomposition. Auder et al. (2010) |fT8l| extended this approach by 
preliminary classification. 

In this work, we propose a new functional approach involving a combination of Shape Invariant Model (SIM) 
approach and Gaussian Process modelling. Considering / time steps and a training set X" = {xi, .., x„} the simulator 
output consist then in a set of curves: 



Y« = [Yij = F(xi,tj), l<i<nA<j<J 

In our practical example (the CO2 storage case study) we have observed that usually these curves have quite 
similar behaviour. So that, we assume that there is a mean pattern from which all the curves can be deduced by a 
proper parametrical transformation. The shape invariant model assumes that the curves have a common shape, which 
is modelled nonparametrically. Then curve- specific diff'erences in amplitude and timing are modelled with the linear 
transformations such as shift and scale |[T9l[2Qll2TI . Hence, we consider the following model: 

Ytj = a; f(tj-o;)^v;^(T;stj, i<i<nA<j<j d) 

where tj stands for observation times and 6* = {6i, 1 < / < n}, v* = {v/, 1 < / < n}, a* = {ai, I < i <n} eW are vectors 
of transformation parameters corresponding to horizontal and vertical shifts and scales of the unknown function /. 

Without loss of generality we assume that the simulator provides the data at some constant time intervals. Time t 
is then equispaced in [0, T] with / time steps. The unknown errors cr*Sij are independent zero-mean random variables 
with variance cr*^. Here, we can also assume that cr*^ = 1 without loss of generality. 

The approach proposed in this work for the functional outputs modeling combines an efficient estimation of the 
transformation parameters (a*, ^*, v*) and a subsequent GP modeling for these parameters that can be used to predict 
new curves without running the simulator. 

The paper is organized as follows. Section 2 presents the basics of GP modeling and the model validation criteria. 
Section 3 describes the method of efl&cient estimation of transformation parameters in the shape invariant model. The 



method is illustrated with an example on an analytical function. In Section 4 we present the forecast algorithm for 
dynamic simulators. Section 5 provides the practical application of the algorithm with a CO2 storage reservoir case. 

2. GP based metamodeling 

The method proposed in this paper is based on a combination of a shape invariant model and a Gaussian process 
metamodel. In this section, we recall the basics of GP modeling or kriging. 

The idea of modeling an unknown function by a stochastic process was introduced in the field of geostatistics by 
Krige in the 1950's |22| and formalized in 1960's by Matheron (1963) 1 10|. Later Sacks et al. (1989) |8| proposed 
the use of kriging for prediction and design of experiments. The theory and the algorithms are formalized in O, 1231 
and 13. 

Consider the output of a simulator as an unknown deterministic function ¥(x) e R, where x g Q c R^ is a 
specified set of selected input parameters. The function F is only known in predetermined design points: X" = 
{(xk, Fk = F(xk)) ,1 <k <n}. The objective is to predict the function Fq = F(xo) for some new arbitrary input xq. The 
function is modeled as a sample path of a stochastic process of the form: 

m 

¥(x) = 2 hj(x) . fij + Z(x) = lf\i{x) + Z(x) (2) 

where: 

• p\i{x) is the mean of the process and corresponds to a linear regression model with preselected given real- 
valued functions h = {hi, 1 <i<m}. Here, we only consider the case h = 1. 

• Z(x) is a centered Gaussian stationary random process. It is defined by its covariance function: C{x,y) = 
E[Z(x)Z(y)] = cr^R(x,y). R(x,y) is the correlation function and cr^ = E[Z(x)^] denotes the process variance. 
Stationarity condition assumes: R(x,j) = R(|x-}^|), where \x-y\ denotes the distance between x e fl and 

Numerous families of correlation functions have been proposed in the literature. We use here Gaussian correlation 
function, the special case of the power exponential family. The power exponential correlation function is of the 
following form: 



f d /v.-.^YM 



R(x,3;) = exp 
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(3) 



where dj = \xj-yjl < pj <2 and 6j > 0. The hyperparameters (^1, .., 0^) stands for correlation lengths which aff'ect 
how far a sample point's influence extends. A high Of means that all points will have a high correlation (F(x/) being 
similar across our sample), while a low Oi means that there are significant diff'erence between the F(x/)'s |24|. The 
parameters pj corresponds to the smoothness parameters. These parameters determine mean-square diff'erentiablity of 
the random process Z(x). For pj = 2 the process is infinitely mean-square diff'erentiable and the correlation function 
is called Gaussian correlation function. Hence, Gaussian correlation function is infinitely mean-square diff'erentiable 
and it leads to a stationary and anisotropic process Z(x) |[9l[23|. Regardless the choice of a correlation function, 
the estimation of hyperparameters (^1, ..,6d) is crucial for reliable prediction. We are using maximum likelihood 
estimation algorithm |9 1 that we will discuss later. 

The experimental design points are selected in order to retrieve most information on the function at the lowest 
computational cost. The number of design points for a reliable response surface model depends on the number of 
inputs and on the complexity of the response to analyze |7, 9|. Latin Hypercube Designs (LHD) provides a uniform 
coverage of the input domain. If we wish to generate a sample of size n, first, we partition the domain of each 
variable in n intervals of equal probability. Then, we randomly sample xi according to the distribution of each of the 
n intervals. Further, for each of the n values for xi , we randomly select one interval to sample for X2 , so that only one 
sample of X2 is taken in each interval. We continue the process of a random sampling without replacement until all the 
variables have been sampled. As a result we generate a sample where each of d inputs is sampled only once in each 



of N intervals. Latin hypercube designs have been applied in many computer experiments since they were proposed 
by MckayetaL, (1979)17 1. 

In this work, we use modified version of LHD - maximin LHD. It is based on maximizing a measure of closeness 
of the points in a design D" : 

max min d(xi,X2) 

design D" xi,X2gD" 

It can guarantee that any two points in the design are not "too close". Hence, the design points are uniformly spread 
over the input domain. 

Consequently, when we have the experimental design X" = (xi , .., x„) and the observation data Y" = (F(xi ), .., F(Xn)) 
the multivariate distribution according to the model ^ for the Gaussian correlation function can be expressed as: 



y. I-M. 



h^(^o) \o ^2! 1 r^xo) 



H r "^ \ r(xo) R 



where R = (R(x/, Xj))i<ij<n g M"^" is the correlation matrix among the observations; r(xo) = (R(xi, xq), .., R(x„, xq))^ g 
R" is the correlation vector between the observations and the prediction point; h^(xo) = \hj(xo)\ g R'^ is the 

vector of regression function at the prediction point xn and H = (hj(xi)) e R^x'^ is the matrix of regression 

\ -^ / l<i<n,l<j<m 

functions at the experimental design. The parameters P and a are unknown. 

Considering the unbiasedness constraint, the parameter P is replaced by the generalized least squares estimate P 

in ([2]). Here, P is of the following form: P = (H^R"^HJ H^R"^Y". Assuming that the correlation function is the 
Gaussian correlation function, the prediction is therefore given by: 

F(xo) = h^(xo) • P + f (xo) R-^ (f" - H . p) (4) 

The hyperparameters 6 = (^1, .., 6ni) and the process variance a^ are estimated by Maximum Likelihood (MLE). 
Using the multivariate normal assumption, the MLE for a^ is: 

0-2 = 1 (y« - HP) R-^ (y" - HP) (5) 

Knowing the estimations forp and &^, the coeflftcients 6 are estimated by maximizing the log likelihood: 

/ ()&, (t^ 6>) = - - [m log &^(6) + log(det(R(e)) + m] (6) 

The function ^ depends only on 0. 

After having estimated all the model parameters, we now need to validate the model. In this work for estimation 
of prediction accuracy of the model, we use the predictivity index, Q2, and root mean squared error, RMS E . The 
predictivity index is calculated basing on cross validation and it has the following form: 



RMSE := 



\/j ^2 := 1 -r- (/) 

^tt ^ I.U{F(xi)-Ff 



Sx/xi i^ the kriging model computed using all the design points X" excepting Xf and F is the mean of F(xi). 

The closer Q2 is to 1 or RMSE is to 0, the higher is the model predictivity. These criteria can be also calculated 
on a separate validation test data by performing additional simulations. It provides higher accuracy measure though it 
requires additional time costs. 

3. Shape Invariant Model 

In this section, we discuss the shape invariant model representation and the procedure for efl&cient parameters 
estimation. 



The shape invariant model was introduced by Lawton et al. (1972) 1 19|. The model assumes that we are working 
with a set of curves that have a common shape function that is modeled nonparametrically. The deformation of this 
function is modeled parametrically by choosing a proper parametrical transformation. We consider a class of linear 
transformations only. These parameters can be normally interpreted as shift in timing or scale in amplitude. For 
this reason, shape invariant model is widely applied to study periodic data such as temperature annual cycle |'25'| or 
traffic data analysis |26|. Indeed, in these cases there is always some variability in time cycles or amplitude. The 
model has also been used to study biological data , where the departure from the pattern can be caused by a different 
environmental conditions 112111201 . In this work, we use the model to propagate uncertainty on a reservoir simulator 
output data. For example, we consider a reservoir pressure during CO2 injection, then shift in time is caused by 
different moment of stopping injection. So that, by applying shape invariant model, we can study the influence of the 
model input parameters on the overall shapes of the selected output. 

In figure ([T]) we display three possible transformations that we consider in our work: horizontal shift 1(a) vertical 
shift|l(b)|and vertical scaling pX^)] The bold line represent original pattern shape. 
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(a) Horizontal Shift (b) Vertical Shift 

Figure 1 : Parametrical transformation examples 



(c) Vertical Scale 



We are interested in the common shape as well as in the efficient estimation of transformation parameters. So that, 
we can reproduce any curve and align it to the pattern. Moreover, we can make a prediction of a possible new curve 
for an input configuration xq by modeling the transformation parameters. 

As already mentioned, for an experimental design X" = {xi, .., x„} we have a set of observations: 

Y" = lYij = F(x/, tj), I < i <n, 1 < 7 < /}, where x/, I < i < n is the set of preselected input parameters and 
tj, I < j < J refers to the time sample. So that, Ytj is the /^ observation on the f^ experimental design unit, 
with I < i < n and I < j < J. Thereby, the idea is to find a general pattern curve that makes possible subsequent 
transformation of any curve to this selected pattern with properly adjusted transformation parameters. 

We focus here on linear transformations. Thus, the model structure may be written as: 



nj = ^rf(tj-o;) + v; + (T''^-st^ 



(8) 



where tj, I < j < J are observation times which are assumed to be known and equispaced in the interval [0, T]. 
The vector of parameters: (a*, 0*, v*) = (ofi, .., Qf„, 61, .., On, vi, .., v„) is unknown as well as the pattern function /. The 
errors Sfj are i.i.d. with a normal distribution, (/, j) e {1, .., ^} x {1, .., /}. It characterizes observation noise. Without 
loss of generality we may assume that cr*^ = 1. The variance does not affect the parameters estimation procedure and 
the method still works for a variance cr*^. The function / is assumed to be 27r-periodic ll25ll26ll . 

In this section we will provide the algorithm for efficient estimation of transformation parameters under unknown 
function pattern /. Since the functional pattern / is unknown, the pattern is replaced by its estimate. So that, it 
seems natural to study the problem ([S]) in a semi-parametric framework: the transformation shifts and scales are the 
parameters to be estimated, while the pattern stands for an unknown nuisance functional parameter. We use an M- 



estimator built on the Fourier series of the data. Under identifiability assumptions it is possible to provide a consistent 
algorithm to estimate (a*, ^*, v*) when / is unknown. The algorithm is described in details in the following subsection 
with an illustration on an analytical function example. 

3.1. Model Assumptions 

Consider Y" = F(X", t) = {F(xi, tj), 1 < i < n,l < j < J} is (n x J) matrix of observations. We model these 
observations in the following way: 



Ytj = a; . f(tj - ^;) + v* + Sij, l<i<nA<j<J 



(9) 



where / : R ^ R is an unknown 27r-periodic continuous function, ^* = (^*,...,^*), a* = (»*,...,»*), v* = 
(vi, ..., v„) G R" are unknown parameters, stj is a Gaussian white noise with variance equal to 1. The time period is 
translated in such a way that: [0, T[^ [0, 2;r[, therefore ti = ^2?:, j = 1, .., / are equispaced in [0, 27r[. 

The objective is to estimate the horizontal shift 0* = (61, ...,^*), the vertical shift v* = (vi, ..., v„) and the scale 
parameter a* = (a^ ..., a*) without knowing of the pattern /. Fourier analysis is well suited for the selected structure 
of the model. Indeed, this transformation is linear and shift invariant. Therefore, by applying a discrete Fourier 
transform to ^ and supposing / is odd, the model becomes: 



^ I ale-^^'lci(f) + Wki, l<k<n, < |/| < (/ - l)/2 
\ alco(f) + V* + Wko, l<k<n, 1 = 



""'' - \ alco(f) + V* + Wko, l<k<n, 1 = ^'""^ 

where c/(/) = j Yjm^i fOm)e~^^^^ , (|/|<(/-i)/2) are the discrete Fourier coefficients and Wki, (i<k<n, |/|<(/-i)/2) is a 
complex white Gaussian noise with independent real and imaginary parts. 

We also notice that in order to ensure the identifiability of the model ( p^ we are working in the parameter space: 
A = {(a%6/%v*) G [-;r,;rp'^: ai = l,6/i = 0,vi = o}. 

To summarize, in this section we estimate the transformation parameters (a*,^*,v*) without prior knowledge of 
the function /. The estimation depends on the unknown functional parameter (c/ (/))|/|<(/-i)/2' the Fourier coefficients 
of the unknown function /. So that, we consider a semi-parametrical method based on an M-estimation procedure. 
M-estimation theory enables to build an estimator defined as a minimiser of a well-tailored empirical criterion that is 
given in the following subsection. 

3.2. Parameters estimation procedure 

The goal is to estimate the vector of parameters (a*, ^*, v*) that depends on the Fourier coefficients of the unknown 
function /. We consider a semi-parametric method based on an M-estimation procedure 1 26 1 . 

To construct an M-estimator, we define the rephased (untranslated and rescaled) coefficients for any vector (^, a, v) g 
A: 

~ . . ^ ( ^^''''*/' ^^k<n, 0<|/|<(/-l)/2 



a^ (dki-Vk), l<k<n, 1 = 



and the mean of these coefficients: 



1 " 
(a, 6,v)=- y Cki(a, 6, v), |/| < (/ - l)/2 
n ^ 



ci 

n 

k^l 



Therefore, for (a*, ^*, v*) we obtain: 

c,/(a%r, V*) = ciif) + ^ e'^'^^wu l<k< 



n 



\^e'^oi-wu 



ci(a*,0%v*) = c,(f)+-y ': 



Hence, \cki(a, 6, v) - ci(a, 6, v)p should be small when (a, 6^ v) is closed to (a*, ^*, v*). 
Now, consider a bounded positive measure yu on [0, T] and set 

Si:= J expJ^JjyuM, (/gZ) (11) 

Obviously, the sequence (Si) is bounded. Without loss of generality we will assume that ^o = and |^/| > 0, / 9^ 0. 
Assume further that 2/ l^/P k/(/)P < +<^- So that f ^ jiis 3. well defined square integrable function: 



f^/d(x) = f f(x - y)d/d(y) 
We construct the following empirical contrast function ([12]): 



J -I 

n 2 



-1 " 2 

M„(a, ^, V) = - . ^ ^ l^/p \cu{a, 6, V) - c/(a, 6, vf (12) 

^ 2 

The random function M„ is non negative. Furthermore, its minimum value is reached closely to the true parameters 
(or*, ^*, V*). We define the estimator by: 

(a,6,'v)n = arg min Mn(a,6,v) 

a,6,veA 
— P 

The proof of convergence (a, 6,'v)n — ^ (a*, ^*, v*) and asymptotic normality of the estimators can be found in 



1261 [251. The weights 61 in \\A) are chosen to guarantee the convergence of the contrast function to a deterministic 
function M„ and to provide the asymptotic normality of the estimators. Moreover, the convergence can be speeded up 
by proper selection of weights. The analysis of covergence at diff'erent weights is presented in |26| . In this work we 
used the weight di = l/\lf withyS = 1.5. 

The computation of the estimators is very fast since only a Fast Fourier algorithm and a minimization algorithm 
of a quadratic functional are needed. The procedure is summarized in Algorithm ([T]) 

Algorithm 1 Parameters estimation procedure 

Input: Input set of curves from experimental design Y" = {YijJ = 1, .., ^; j = 1, .., /} 

Output: Transformation parameters estimation (a*, ^*, v*) 

Define the identifiability condition: A = {(a%6/%v*) e [-;r,;r[^^'^: ofi = l,^i = 0, vi = o} 

Compute the matrix of discrete Fourier coefficients D = [dku k = 1, .., ^; |/| < (J - l)/2} 

Compute the matrix of rephased Fourier coefficient C = {cki, k = 1, .., ^; |/| < (J - l)/2} 

Compute the vector of mean of rephased coefficients C = {c/, |/| < (J - l)/2} 

Choose the weight sequence Si 

Define M,(a, 0,v) = '- - ZLi T^T i^ l^/l' \^ki(a^ ^' ^) " ^/(^' ^' ^f 

Compute (a,^, v) = arg mino^^^^veA M„(a,^, v) e R^x^^-^) 

Return: (a,?, v) g R3><("-i> 



3.3. Analytical Function Example 

In this section the shape invariant model and the efficient parameters estimation are presented on an analytical 
function. The minimization algorithm used in the estimation procedure is a Newton-type algorithm. 
Let us consider the following function: 



/(x) = 20.(1 -x/(27r)).x/(27r) 



Simulated data are generated as follows: 

Yij = a]- f(tj -6;) + v;+stj, l<i<nA<j<J 

with the following choice of parameters: J = 5, N = lOl, tj = ^2?:, 1 < 7 < / are equally spaced points on 
[0, 27r[. Transformation parameters (^*, a*, v*) are uniformly distributed on ]0, 1], where ^^ = 0, v^ = 0, al = I; the 
noise Sji,j= 1, ..,/,/ = 1, .., A/^ are simulated with a Gaussian law with mean and variance 0.5. 

Results are displayed in Figure[2] The function / is plotted by a solid red line in Figure 2(d) Figure 2(a) [shows the 



original simulated noisy data Yij. The cross-sectional mean function of these data is presented in Figure 2(d) by black 



dotted line. Figure [2(b)] plots estimated transformation parameter versus the originally simulated parameters. As it 
can be seen, the estimations are very close to the original parameters. The inverse transformation using the estimated 
parameters is displayed in Figure [2(c)] and the mean function of restored curves is displayed in Figure [2(d)] by blue 



dashed line. Figure [2(d)[ compares the cross-sectional mean of inversely transformed data and the cross-sectional 
mean of originally simulated data. Despite the noise, it is noticeable that the data after inverse transformation are 
much more closer to the original function / than the original cross-sectional mean function. 

This analytical example shows that the method is effective in estimating the transformation parameters of the 
shape invariant model. In the next section, we will explain how this model can be applied in reservoir engineering 
forecast problems. 

4. Functional data approximation 

To apply the shape invariant model approach to approximate functional data from a dynamic simulator, firstly 
we have to modify the parameters estimation procedure for large number of curves. When we are working with 
uncertainty modeling, we always start from an experimental design X" and a set of observation Y". As we have 
mentioned, the number of design points depends on the number of inputs and on the complexity of the response. 
So that, some optimization problems could arise when we compute the contrast function with a large number of 
curves. It can be time consuming and the results can be inaccurate. Therefore, we propose the following mod- 
ification to the Algorithm ([Tj): the original observation data is split into blocks and the optimization procedure is 
then performed on each of the blocks. Following the identifiability condition, we are working on a compact set 
A = |(a*, 6*, V*) e [-n, 7r[^^": ai = l,6i = 0, Vi = Ol. This condition should be satisfied on every block optimization 
by adding the first reference curve to the block. 

Algorithm 2 Parameters estimation procedure for large n 

Input: Input set of curves from experimental design Y" = {YijJ = 1, .., ^; j = 1, .., /} 
Output: Transformation parameters estimation (a*, ^*, v*) 
Split the observation data into N^ blocks of (K -\- 1) curves 
form = l,..,Nb do 

Define block curves Y^^^ = {7i, Y^m-i)iK^i)...,Y„J = [Yij, i = 1, .., /^ + 1, j = 1, .., /} 

Perform Algorithm ^ 

Compute (a, ^, v) = arg min M„(a, 6, v), where (a, O^v) e R^^^ 

a,6,veA 

end for 

Return: (a,?, v) g R^^^^^-^^ 

With this procedure, we do not have limitations on experimental design of any size. As soon as we have estimated 
the parameters for every curve from observation data set, we can formulate the prediction algorithm. Instead of 
reproducing the simulator output for a prediction point xq at every time step, we model the whole output curve 
with the transformation parameter. This curve will provide the approximation of the output for the selected input 
configuration xq at each of considered time steps {tj, 7 = 1, .., /}. The transformation parameters for the input xq are 
evaluated with the Gaussian process response surface modeling. The model is based on the experimental design and 



the set of evaluated transformation parameters calculated for the observation data curves. The prediction framework 
for an arbitrary input configuration xq is presented by the following Algorithm ([3]). 

Algorithm 3 Prediction algorithm for dynamic simulator 

Input: Dynamic simulator Y = F(x, t) with t ^{tj, j = 1, .., /} and prediction point xq 

Output: Prediction Y^ = F(Jco, tj) for all j = 1, .., / 

Generate an experimental design X'^ = (xi, .., x„) to span the space of interest 

Evaluate Y" = F(X", tj) at every time step tj, I < j <J 

Generate a set of discrete curves {Ytj}, / = 1, .., ^; j = 1, .., / 

Estimate the (a, 6, v) g (Wf with Algorithm[2] 

Construct new experimental design for the function of parameters: (X", ^(X")), (X", a(X")) and (X", v(X")) 

Estimation of hyperparameters for GP models of transformation parameters 

a(xo), 0(xo) and v(xo) are approximated with corresponding GP models 

Reproduce: F(xo, tj) = a(xo)f(tj - 6(xo)) + v(xo) for all {tj, j = 1, .., /} 

Return: Discrete time series Y^ = F(xo, t) with t e {tj, j=l,..,J} 

Summing up, with this proposed algorithm the problem of response surface modeling for dynamic simulators is 
reduced from single step GP modeling for each of / time steps to an optimization problem and a GP modeling for 
the transformation parameters. However, it is worth to mention that before performing the algorithm it is important to 
analyze the curves behaviour for the observation data set. Studying the curves characterization, probably we may fix 
for example vertical shifts v or horizontal shifts 6 at zero. Also, if the curves have significantly diff'erent behaviour at 
diff'erent time intervals we can split the observation data in time as well in order to achieve higher prediction accuracy. 

The next section presents an application with a dynamic reservoir simulator case. 

5. CO2 storage reservoir case 

Carbon Capture and Storage technology stands for the collection of CO2 from industrial sources and its injection 
underground. Carbon dioxide is stored in a deep geological formation that is sealed on a top by a low permeability 
cap rock formation. Subsurface storage of CO2 is always associated with an excess pressure in the reservoir and one 
of the primary environmental risks is a pressure-driven leakage of CO2 from the storage formation. In order to assess 
the risk of CO2 leakage through the cap rock we consider a synthetic reservoir model. The model is made up of three 
zones ([3]): 

• a reservoir made of 10 layers 

• a cap-rock made up of 9 layers 

• a zone- to- surface composed of 1 layer 



1 layer 



9 layers 



10 layers 




PHI = 0.25 
K = 200md 



PHI = 0.37 
K = 10-s md 



PHI = 0.25 
K = 200md 



Figure 3: Reservoir Model 

The XY size of the grid is set at 10 km total length. Each layer is 5m thick, including the cell above the cap-rock. 
The total number of cells is 13520 (26x26x20 model grid). The structure of the reservoir is reduced to its simplest 
expression. The zone above the cap-rock (up to the surface) is currently set to 1 layer. The salinity of the water is 
35gm/l. The temperature of the reservoir is set to 60C and the initial pressure is hydrostatic. The injection bottom rate 
is set to 10E-h06 tons/year. The fracture pressure is estimated by geomechanical experts to 158 bars. Exceeding this 
value of reservoir pressure can lead to a leakage. The simulation period is 55 years that include an injection period 
of 15 years followed by 40 years of storage. In this study we analyze the possibility of leakage through a cap rock. 
Therefore, we consider pressure in the storage reservoir as an objective function to be approximated. 

The uncertain parameters selected for this study characterize the reservoir and the fluid properties. It implies 
difl'erent CO2 flowing possibilities between the reservoir layers. The distribution law for the parameters is uniform. 
Table ([T]) represents the parameters description with their range of minimum and maximum values. 



Name 


Description 


Min Max 


PORO 

KSAND 

KRSAND 


Reservoir Porosity 

Reservoir Permeability 

Water relative permeability end-point 


0.15 0.35 
10 300 
0.5 1.0 



Table 1 : Uncertain Parameters 



We start from the observation data Y" - [30; 55] matrix of simulator outputs. By means of Algorithm ^ we 
provide the transformations parameters estimations. Figure ( |4(a)| ) provides the original set of curves and Figure ( |4(b)| ) 
represents the same set after inverse transformation. The pattern curve is diflferentiable after the inverse transformation 
with the estimated parameters. 
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Reservoir Pressure Development 



After Inverse Transformation 





30 
Years 



(a) Observation data (b) Restored data 

Figure 4: Original observation data and data after inverse transformation 

As we are sure here that the model parameters are efficiently estimated, we can proceed with the next step of 
prediction algorithm ([3]). The next step is to build Gaussian process response surface models for the transformation 
parameters basing on the estimations. 

In figure ([5]) we display the model validation criteria: Root Mean Square Error (RMSE) and predictivity indices 
(Q2), calculated separately for every year. The criteria were computed with the help of additional confirmation test 
data. The low predictivity in the first and last years is caused by low variance of data in that period. In general, the 
method provides reliable level of predictivity. It is also reflected by crossplots of test and predicted data. Figure ( [6(a)| ) 
is based on new proposed approach and Figure ( |6(b)| ) corresponds to single step GP modelling. Both methods provide 
a good level of approximation. 



11 



Crossplot of Predicted and Test Data, SIM method 



Crossplot of Predicted and Test Data. Single step GP 




110 115 

SIM predicted values, bars 

(a) SIM method 




110 115 

SIM predicted values, bars 

(b) Single Step GP 



Figure 6: Crossplot comparison 

Table ([2]) compares the simulation CPU time for both methods. SIM method provides suffucuently accurate results 
(although less accurate than the single step GP approach), but with a CPU reduction of a factor five. 





SIM method 


Single step GP 




Optimization 


Parameters modeling 


Total 


CPU time 


00:00:35 


00:00:15 


00:00:50 


00:04:28 



Table 2: CPU time comparison 

It is worth to mention, that in this study we consider a simple model with only 3 uncertain parameters. So that, to 
estimate the function with GP model at every single step takes approximately 10 seconds. For more complex functions 
and more input uncertain parameters involved a single step model evaluation can take up to 10-20 minutes. So that, 
for the same simulation period of 55 years CPU time can increase to 10-20 hours. Whereas SIM approach does not 
depend on number of time steps and the method always conclude only a single optimization problem and as maximum 
3 GP models for the transformation parameters. 

6. Conclusion 

This paper focuses on two general problems. First, we introduce a shape invariant model approach and we provide 
an efficient algorithm for estimation of transformation parameters of the model with the method specification for 
large sets of curves. Second, we suggest the application of this approach to model the time series outputs from a 
dynamic simulator. The proposed method reduces the problem of functional outputs modeling to one optimization 
problem and three GP response surface models. We have tested the method with a CO2 storage reservoir case. We 
have also compared the method with the standard single-step approach. Presented numerical results show that the 
method provides satisfactory and comparable predictivity at lesser CPU time. The method also does not depend on 
the number of the involved time-steps. It can be very advantageous when we are working with a model involving large 
number of times steps such as CO2 storage when the reservoir model simulation period can include up to hundreds or 
thousands time steps. However, if the set of output curves have significantly diff'erent behaviour, preliminary curves 
classification may be required. 
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Figure 2: Analytical Example 
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Root Mean Squared Error 
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Figure 5: Predictivity Indices and RMSE 
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